﻿#include <QCoreApplication>

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include "vector.h"
#include "convolution.h"
#include "iir.h"

using namespace std;
using namespace splab;

// 实现零相位滤波函数
template<typename Type>
Vector<Type> filtfilt(const Vector<Type>& signal, const Vector<Type>& numCoefs, const Vector<Type>& denCoefs) {
    // 正向滤波
    Vector<Type> forward = convolution(signal, numCoefs);
    forward = convolution(forward, denCoefs);

    // 反转信号
    Vector<Type> reversed = reverse(forward);

    // 反向滤波
    Vector<Type> backward = convolution(reversed, numCoefs);
    backward = convolution(backward, denCoefs);

    // 再次反转得到最终结果
    return reverse(backward);
}

// 生成测试信号：正弦波加高频噪声
template<typename Type>
Vector<Type> generateTestSignal(int length, Type fs) {
    Vector<Type> signal(length);
    Type t;

    for(int i = 0; i < length; i++) {
        t = Type(i) / fs;
        // 生成1Hz的正弦波加上10Hz的高频噪声
        signal[i] = sin(2 * PI * t) + 0.5 * sin(2 * PI * 10 * t);
    }

    return signal;
}

int main(int argc, char *argv[])
{
    QCoreApplication a(argc, argv);

    // 设置采样频率和信号长度
    double fs = 100.0;  // 采样频率100Hz
    int signalLength = 1000;  // 信号长度1000点

    // 生成测试信号
    Vector<double> signal = generateTestSignal<double>(signalLength, fs);

    // 设计低通IIR滤波器（Butterworth）
    IIR filter("lowpass", "Butterworth");
    filter.setParams(fs, 2.0, -3.0, 5.0, -40.0);  // 截止频率2Hz
    filter.design();

    // 获取滤波器系数
    Vector<double> numCoefs = filter.getNumCoefs();  // 分子系数
    Vector<double> denCoefs = filter.getDenCoefs();  // 分母系数

    // 应用零相位滤波
    Vector<double> filteredSignal = filtfilt(signal, numCoefs, denCoefs);

    // 输出结果
    cout << "原始信号和滤波后信号的前20个采样点：" << endl;
    cout << setw(10) << "时间(s)" << setw(15) << "原始信号" << setw(15) << "滤波后信号" << endl;
    cout << "------------------------------------------------" << endl;

    for(int i = 0; i < 20; i++) {
        cout << setw(10) << fixed << setprecision(3) << i/fs
             << setw(15) << signal[i]
             << setw(15) << filteredSignal[i] << endl;
    }

    return 0;

    return a.exec();
}
